Phase Stability in the Two dimensional Anisotropic Boson Hubbard Hamiltonian 
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The two dimensional square lattice hard-core boson Hubbard model with near neighbor 
interactions has a 'checkerboard' charge density wave insulating phase at half-filling and sufficiently 
large intersite repulsion. When doped, rather than forming a supersolid phase in which long range 
charge density wave correlations coexist with a condensation of superfluid defects, the system instead 
phase separates. However, it is known that there are other lattice geometries and interaction patterns 
for which such coexistence takes place. In this paper we explore the possibility that anisotropic 
hopping or anisotropic near neighbor repulsion might similarly stabilize the square lattice supersolid. 
By considering the charge density wave structure factor and superfluid density for different ratios of 
interaction strength and hybridization in the x and y directions, we conclude that phase separation 
still occurs. 

PACS numbers: 71.10.Fd, 02.70.Uu 



I. INTRODUCTION 

In fcrmionic systems the competition between 
phase separation and (non-phonon mediated) 
superconductivity has been a long-standing and 
active area of research. On the one hand, Cooper pair 
formation requires an effective attractive interaction, but 
on the other hand this attraction might 'go overboard' 
and lead to the agglomeration of larger clusters of 
particles. This latter possibility is especially likely in 
commonly studied models of superconductivity like the 
Hubbard and t-J Hamiltoniansi~— which do not include 
a long-range Coulomb interaction that opposes phase 
separation. 

A very similar set of issues arises in the interplay 
of Bose-Einstein condensation (and the closely related 
phenomenon of superfluidity) with charge order in 
bosonic systems, where the possibility of coexistence of 
diagonal and off-diagonal long range order is termed a 
'supersolid' (SS) phased. Bosonic systems are commonly 
studied in the continuuin&i -, motivated by liquid 4 He. 
However, as with fermions, issues such as supersolidity 
can also be addressed within a tight-binding lattice model 
governed by the boson Hubbard Hamiltoniaraii, 
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Here cj, a i , and n i are boson creation, destruction, and 
number operators. Uj is an intersite hopping between 
sites i and j. In this paper we consider a square lattice 
of linear size L and total number of sites N — L 2 . Ui is an 
on-site repulsion, and Vij is an intersite repulsion. The 
chemical potential controls the filling of the lattice. 



Most studies of the boson-Hubbard Hamiltonian have 
considered the case of translationally invariant and short 
range interactions Ui = U, Vij = V\ (V2) for near neighbor 
(next near neighbor) pairs ij, and near neighbor hopping 
tij = t. More complex choices for the Hamiltonian 
parameters are increasingly being studied, most notably 
situations in which the chemical potential models a 
quadratic trap to confine the bosons 1 ^—. 

In the uniform interaction, hard-core limit, U = 00, a 
site can be either empty or singly occupied. When a near 
neighbor, Vi, is turned on at half-filling, the bosons begin 
to occupy only one of the two sublatticcs into which the 
square lattice is divided, owing to its bipartite character. 
For V\ sufficiently large and at low temperature, this 
tendency becomes sufficiently pronounced to result in 
long range order: the boson occupations alternate in 
a checkerboard pattern, with a charge ordering vector 
q = (n, 7r) and charge correlations extend to large 
separation. Doping away from half-filling creates defects 
in this charge density wave (CDW) phase. If the defect 
concentration is sufficiently low, it is possible that long 
range charge order is not destroyed, even though the 
defects can move through the ordered lattice and form 
a superfluid. In this way a supersolid would be realized. 

However, it was show n 16 ' 17 that this coexistence 
does not occur for hardcore bosons on the 2D square 
lattice with near neighbor interactions. Rather, there 
is a thermodynamic instability to phase separation in 
which different spatial regions are superfluid and charge 
ordered. Despite the absence of simultaneous diagonal 
and off-diagonal long range order in this most simple 
scenario, other lattice geometries have been shown to 
exhibit supersolid phases, and, indeed, by now there is a 
considerable numerical literature as a function of lattice 
geometry, dimensionality, filling and interaction 1 ^ - — . 
More exotic processes within the Hamiltonian, such as 
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ring exchange, are also known to play a potential role in 
supersolidit y 38 ' 39 . 

One of the lines of investigation into supersolidity 
has pursued the effect of anisotropy. Thus, although 
the checkerboard supersolid with q = (tt, tt) ordering 
wavevector is unstable, it was also shown that a striped 
supersolid is stable^ 7 -. Such a phase is caused by a next- 
near neighbor interaction V2 which drives charge order 
at wavevector q = (tt,0) or (0,7r) which is then doped 
to allow for mobile defects. In this case the underlying 
Hamiltonian is isotropic, but a spontaneous breaking of 
the Z 2 symmetry associated with order along the x or y 
directions occurs. Recent studies have also shown that 
symmetry breaking in the Hamiltonian itself can give rise 
to striped supcrsolids. In particular, Chan et al^ argued 
that an attractive interaction in one spatial direction and 
a repulsive interaction in the other will cause supersolid 
formation. 

These known instances of stability of the striped 
supersolid, as well as the general delicacy of the energy 
balance involved in whether phase separation occurs or 
not£, suggest that a checkerboard supersolid might be 
stabilized by anisotropic hoppings t x ^ t y or interactions 
V x / V y in the x and y directions. This is the topic of 
the present paper. Stable supersolid phases have already 
been found for soft-core bosons in one, two and three 
dimension a 22 ' 23 i 36 ' 37 i 41 . When the contact repulsion, U, 
is large compared to the near neighbor interaction, V±, 
(the mean field value is U > 4Vi in two dimensions), 
the system behaves as in the hard core case and phase 
separates when doped away from half filling. When V\ 
dominates, U < 4Vi, and the system is doped above 
half filling, the extra bosons will go on already occupied 
sites. These doubly occupied defects then delocalize 
and yield a supersolid phase. On the other hand, if 
the system is doped below half-filling, it will typically 
undergo phase separation. Remarkably, however, there 
is a very narrow parameter range where the system will 
become supersolid when doped below half fillin g 36 ! 37 . 
This history of successful searches for supersolidity in 
soft-core models prompts us to focus here on the hard- 
core case. 

After introducing our Quantum Monte Carlo 
methodology and the observables used to characterize 
the phases, we show results indicating that supersolids 
do not emerge from a checkerboard charge ordered 
pattern driven by near-neighbor interactions, even when 
those interactions or the hoppings are made anisotropic. 
The dependences of the supcrfluid density p s and the 
CDW structure factor S(q) on V\ x jV\ y and t x /t y are 
calculated. 



works in continuous imaginary time (avoiding the need 
to extrapolate to zero discretization mesh of the inverse 
temperature /3) . Here we use a recent new formulation of 
the method which includes global space-time updates^ 
to explore phase space efficiently and to sample the 
winding and hence supcrfluid density. We use a canonical 
formulation which allows us to work with fixed particle 
number, although the method can also be implemented 
in the grand canonical ensemble. Note that even within 
a canonical implementation, it is still possible to obtain 
the chemical potential as a finite difference of ground 
state energies at adjacent particle numbers, p(n) = 
E (n + l)-E (n). 

We will characterize the phases of the Hamiltonian 
Eq. [1] by examining the energy, E = (H) and also the 
supcrfluid density p s and CDW structure factor S(q). 
The supcrfluid density is obtained using the procedure 
described by Pollock and Ceperley^, in which p s is 
expressed in terms of the winding W of bosonic world 
lines across a lattice edge. 

„, _ (52 ,„ . <2S> ( 2 ) 

Psx 2(3 t x l sy 2/3 t y { ' 

Here we have allowed for the possibilty of inequivalent 
spatial directions, and measure the x and y superfluid 
densities separately. 

The CDW structure factor is the Fourier transform of 
the density-density correlation function. 

s («) = ]^E e< *'Wi+«> ( 3 ) 
3 1 

Here N is the number of sites. In order to establish long 
range order, both p s and S(q) must be extrapolated to 
the N — > 00 limit via finite size scaling. 

The choice of how to compare anisotropic lattices with 
isotropic ones is somewhat arbitrary. Here wc have 
chosen to keep t x + t y fixed as we vary t x away from 
ty, rather than, for example, fixing t y and reducing t x 
to zero. Our rationale is that this preserves the overall 
noninteracting bandwidth W = A(t x + ty). Similarly, 
when we vary the interactions we do so by fixing the sum 
Vix + Vi y . The motivation behind this choice is that in 
the strong coupling limit the energy cost to add a particle 
to a perfect checkerboard solid, 2(Vi x + Vi y ) is preserved. 
That is, the "base of the insulating lobes"— at t = is 
the same for all degrees of anisotropy. 

III. RESULTS 



II. METHODOLOGY 

Our computational approach is the stochastic Green 

function (SGF) approachi 2 ^— . SGF is a finite 

temperature Quantum Monte Carlo algorithm which 



A. Interaction Anisotropy 

We begin by considering t x = t y but V\ x 7^ V\ y . 
Figure Q] shows the checkerboard structure factor S(ir, tt) 
and superfluid density p s as functions of filling p. In 
Figure HJa), S(tt,tt) is robust to anisotropy, declining 
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FIG. 1: (Color online) Checkerboard structure factor (a) 
and superfluid density (b) as functions of filling, comparing 
anisotropic near neighbor interaction strengths Vi x , Vi y with 
isotropic ones. S(tt, 7r) is relatively unaffected by shifting 
Vix away from Vi y : the data for Vi x = 1, Vi v = 5 declines 
somewhat more rapidly with doping away from p = 1 but 
looks qualitatively quite similar. Only when Vi x = does 
the structure factor change behavior significantly. Similar 
remarks apply to the superfluid density. Superflow p 3y in the 
y direction is larger than that in the x direction. Most data 
are for 8x8 lattices at j3 = 16, but lines without symbols are 
for 1 6xf 6 lattices at ft — 32 and show that finite spatial lattice 
and temperature effects are minimal. Here and in subsequent 
figures error bars are smaller than the symbol size and hence 
are not shown. 



only modestly from its value at V\ x = 3, Vi y = 3 to 
its value at V\ x = l,V\ y = 5. However, when V\ x = 0, 
S(ir, 7r) is greatly altered. In this limit, if the hopping 
t is neglected, the lattice decomposes into independent 
chains each with density order at q y = it but arbitrary 
q x . Figure [1] indicates that nonzero V\ x locks these chains 
into place spatially (selecting out q x = tt as well). This 
is likely due to the enhanced kinetic energy when hard 
core bosons in one chain find empty sites on neighboring 
chains. This robustness of checkerboard order in the 
V\ x = limit is confined to p = 0.5. The structure factor 
falls off much more abruptly with doping than when V\ x 
remains nonzero. 

The behavior of p s in Figure [Ub) is also largely 
unaffected by anisotropy. The largest change is an 
enhancement of p s in the vicinity of the insulating solid 
phase at half-filling. Surprisingly, the details of the 
choice of interaction strength become unimportant at low 
density (and high density due to particle-hole symmetry). 
The x and y components of the superfluid response, p sx 
and p sy , have the same qualitative dependence on density 
p, although for the largest anisotropy V\ x — 0, p sy is 
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FIG. 2: (Color online) (a) S(ir, tt) vs V lx /V ly at half-filling 
p — 0.5 and lightly doped p = 15/32. The sum Vi x + Vi y = 6 
is fixed. The flatness of the data except in the vicinity 
of Vi x /Viy — emphasizes the minor effects of anisotropic 
interations, especially at half-filling. (b) p 3X and p sy vs 
Vix/Viy at the same two densities. The system is insulating 
for p = 0.5 except strictly in the Vi x = limit. In the 
doped case, p sx shows little change with anisotropy, while 
p sy increases sharply near Vi x — 0. Spatial lattice size is 8x8 
and inverse temperature /3 = 16. 



roughly twice p sx for all densities. A precise picture of 
the interplay of the anisotropics in charge correlations 
and superfluidity is not clear. Since V\ x << V± y it 
is more likely that x neighbors will both be occupied 
(or empty) than y neighbors. Strings of adjacent 
occupied sites of hard-core bosons would seemingly block 
motion perpendicular to them (and hence suppress p sy ). 
Meanwhile, strings of adjacent empty sites would offer 
channels for an enhanced p sy . A complicating feature of 
any attempt to build such a 'blocking/channel' picture 
is the fact that while the hard-core nature of the bosons 
completely precludes certain types of motion, the strong 
Vi y interactions will also profoundly affect the pattern of 
bosonic flow. The small non-zero value of p sy at V\ x = 
and half-filling is a finite lattice effect. (Sec Fig. |4j) 

The detailed evolution of the CDW structure factor 
and superfluid densities with interaction anisotropy 
Vi x /Viy is given in Fig. O Panel (a) emphasizes that 
S(ir, tt) is generally robust to anisotropy at half-filling 
p = 0.5, but can change dramatically for sufficiently 
small Vi x /Viy when the system is (slightly) doped to 
p = 15/32. Panel (b) indicates similarly that p s is much 
more affected by anisotropy in the presence of vacancies. 
Bosons flow more freely in the direction perpendicular to 
that for which the repulsive interaction is smaller. This 
probably represents the low energy cost of lines of bosons 
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FIG. 3: (Color online)The temperature dependence of 
(a) structure factor S(tt,-k) and (b) superfluid density p s 
(averaged over the x and y directions) are shown. We see that 
for ft > 10, these quantities take their ground state values for 
the parameters we are using. 



which have strong charge ordering in the y direction 
slipping past each other due to the minimal energy cost in 
V\x ■ Of course such a picture is a bit naive in the present 
case when the system undergoes phase separation and 
consists of a mixture of superfluid and CDW regions. 

Figure [Ua,b) contained some comparison of data on 
different lattice sizes. Fig. [3] and Fig. @] show the effects 
of temperature T = 1//3 and L more systematically. 
Fig. [3] establishes that the system is in the ground state 
when j3 ~ 10, validating the choice of j3 = 16 used 
to access the ground state in Fig. [TJ The approach to 
the low temperature limit is somewhat more gradual in 
the extreme anisotropic case V\ x = 0. Fig. [4] allows an 
extrapolation to the thermodynamic limit 1/L — > 0. As 
at half-filling, for a doping p = 15/32 the CDW structure 
factor S(tt, it) is nonzero at large L. Further doping 
to p = 13/32 destroys long range density order. The 
anisotropic system is more sensitive to doping than when 
V lx = V ly . 

Figures Q] [21 HI establish the possibility of a supcrsolid 
phase: Non-zero p s coexists with large S(tv, tt) in the 
vicinity of half- filling. However, it is necessary to 
distinguish a true supersolid phase, in which these orders 
coexist as a single stable thermodynamic phase, from 
a lattice in which both orders are present at separate 
spatial locations. This phase separation is most easily 
detected by looking for an anomalous negative curvature, 
d 2 E/dp 2 < 0, in the thermodynamics. Figure [5] shows 
that such negative curvature is present just below p = 
0.5, unless V\ x = 0. 

The jump in slope p = dE /dp at p = 0.5 is the charge 



FIG. 4: (Color online)The lattice size dependences of (a) 
structure factor S{n, n) and (b) superfluid density p s are 
shown. Solid circles are the isotropic case Vi x = Vi y = 3. The 
system retains checkerboard CDW order when the doping is 
small (p = 15/32) but not for p = 13/32. Open squares are 
the extreme anisotropic case V\ x — 0. In that case there is 
order only at half-filling. p s is larger when the interactions 
are anisotropic, but the finite size scaling indicates that it still 
vanishes at half-filling in the thermodynamic limit even when 
V lx = 0. 
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FIG. 5: (Color online) Energy E as a function of 
density p exhibits the characteristic signal of phase separation 
d 2 E/dp 2 < as long as Vi x remains non-zero. The kinks at 
half-filling imply an abrupt jump in chemical potential, that 
is, the presence of a non-zero charge gap. 
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FIG. 6: (Color online) Phase separation is indicated by a 
chemical potential p which decreases as p approaches p = 0.5 
from below. This phenomenon, previously established in the 
isotropic case, is still present with anisotropy, except in the 
extreme limiting case Vi x = 0. Results for 8x8 lattices at 
/3 = 16 and 16x16 lattices at /3 = 32 show finite-size and 
finite temperature effects are negligible both for the size of 
the Mott gap and the negative curvature (phase separation). 



FIG. 7: (Color online) Kinetic, potential, and total energies 
as functions of the filling p. The kinetic energy is minimized 
(in absolute value) at half-filling, as expected for an insulating 
solid. Phase separation appears to be driven primarily by the 
density dependence of the potential energy, whose decrease 
from a maximum at p w 0.43 results in the negative curvature 
in the total energy. 



gap A c . Figure |5] shows the density p function of 
chemical potential p. The kink at half- filling in Fig. [5] 
is reflected in a plateau in the density. As commented 
earlier, fixing V\ x + V\ y makes the length of this plateau 
A c = 2(Vi x + Vi y ) the same for all degrees of anisotropy 
in the strong coupling (t — > 0) limit. Figure [6] indicates 
that finite hopping t — 1 reduces A c from the t = value 
A c = 12 to A c ~ 8. There is only a relatively small 
reduction in A c going from V\ x = Vi y = 3 to V\ x 
1 . Vi y = 5. Setting Vi x = 0, however, reduces the gap to 
A c ~4. 

In Fig. |5] the signature of phase separation is made 
more evident, in the form of a double valued dependence 
of density on chemical potential p = dE/dp. In other 
words, the chemical potential p, is not a monotonically 
increasing function of density p. The usual Maxwell 
construction can be used to ascertain the relative sizes 
of the supcrfluid and solid regions in this situation. The 
instability of the checkerboard supersolid in the isotropic 
case is knownil. Figures [5] and [5] establish that this 
instability persists in the anisotropic case, except in the 
extreme case of vanishing interchain interaction, Vi x = 0. 
However, as seen in Fig. [TJ the structure factor falls 
off rapidly with doping in this situation, so that long 
range density (solid) order is no longer present. Thus the 
absence of phase separation does not imply the existence 
of a supersolid. 

Figure [7] shows the individual components of the 



energy. At p = 0.5 the system is in an insulating 
checkerboard solid phase, and has a corresponding 
minimum in the absolute value of the kinetic energy. The 
potential energy at first increases from the low density 
limit, reflecting the increased likelihood of adjacently 
occupied sites. At p w 0.43, however, the potential 
energy decreases, despite a rising number of particles. 
Presumably this reflects the onset of density order: a 
coherent pattern of occupation in which adjacent sites 
are rarely occupied. Evidently, the decrease of {niUj) 
for near neighbor sites ij associated with solid formation 
overwhelms the 'trivial' increase {riirij) cx (n.;) (rij) 
expected for noninteracting bosons. 



B. Hopping Anisotropy 

Having explored the possibility of supersolids caused 
by anisotropy in the near neighbor interaction we turn 
to the analogous question when t x ^ t y . Figure [5] shows 
the checkerboard structure factor <S'(7r, 7r) and supcrfluid 
densities p sx and p sy in the x and y directions. The 
structure factor is peaked at half-filling but remains 
relatively large in some interval surrounding that density. 
In general p sx > p sy , the superflow in the x direction 
shows a dip at fillings p = 3/8 and p = 5/8. This 
dip becomes more pronounced as the hopping anisotropy 
increases. We will address this issue, which requires some 
challenging QMC simulations and which is not directly 
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FIG. 8: (Color online) Checkerboard structure factor S(n, tt) 
(top), along with superfluid densities p sx and p ay in the x 
and y directions (bottom). S(tv,tt) is large in the vicinity of 
half-filling and decreases with doping. It is almost unchanged 
by hopping anisotropy. The superfluid density vanishes at 
half-filling, and rises with doping. Beginning with p ax — p sy 
for t x = t y , the superfluid density in the x direction shifts 
upward, and that in the y direction downward, as t x /t y 
decreases. Results for 16x16 lattices at /3 = 32 show that 
finite spatial lattice and temperature effects are minimal. 
For the largest hopping anisotropy, p ax has a dip at fillings 
p = 3/8, 5/8. This will be discussed in a subsequent work. 



related to the questions being addressed here^i, in a later 
paper. 

The thermodynamics, shown in E(p) (Fig. [5]) and p{p) 
(Fig. fTUj) show some features similar to those seen in the 
anisotropic interaction case. E(p) exhibits a kink at p = 
0.5 indicating a jump in the chemical potential and hence 
a charge gap, and also negative curvature near p = 0.5 
suggesting phase separation. p{p) demonstrates these 
features even more directly. The gap (plateau in p(p)) is 
unchanged even for a case when the ratio of hoppings is 
altered by a factor of three from t x /t y = 1 to t x /t y = 1/3. 
This is in contrast to interaction anisotropy to which the 
gap is more sensitive (Figure [5]). 

Figure [TT] provides quantitative detail on the evolution 
of the CDW structure factor and superfluid densities 
with hopping anisotropy. For both densities, the data 
for S(tt,tt) are quite flat over the range shown, 1/3 < 
tx/ty < 1. The superfluid densities, especially in the x 
direction show significant quantitative evolution^, but 
qualitatively they are unchanged, indicating superfluid 
behavior for density p — 15/32 and insulating behavior 
for half-filling p= 1/2. 



FIG. 9: (Color online) Energy E(p) for three different 
hopping anisotropies. All have the kink which signals the 
Mott gap, and also negative curvature (phase separation) in 
the vicinity of half-filling. 
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FIG. 10: (Color online) Dependence of filling on chemical 
potential p(p), extracted from energy data of Fig. [5] At half- 
filling, a plateau in p(p) signals the charge-ordered insulator, 
with adjacent densities showing a thermodynamic instability 
to phase separation. Data for 16x16 lattices at ft = 32 
show that finite spatial lattice and temperature effects are 
minimal both for the size of the Mott gap and the negative 
curvature. The charge gap and negative curvature are largely 
independent of the degree of hopping anisotropy, even though 
ti y /t\ x = 3 is fairly large. Some of this invariance is due 
to our choice t\ x + ti y = 2, which keeps the noninteracting 
bandwidth fixed. 
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FIG. 11: (Color online) (a): S(ir, it) vs t x /t v at p = 0.5 
and p = 15/32. The structure factor is even less sensitive to 
hopping anisotropy than to interaction anisotropy, as might 
perhaps be expected since charge order is directly driven by 
V as opposed to t. (b): p sx and p sy vs t x /t y at p — 0.5 and 
p = 15/32. The sum t x + t v = 2 for all data here. 



IV. CONCLUSIONS 

The simultaneous presence of diagonal and off-diagonal 
long range order, as opposed to a transition between 
distinct ordered phases, is by now well-established in 
many bosonic models. In this paper we have shown 



that this coexistence is thermodynamically unstable in 
a broad range of square lattice hard-core Hamiltonians 
with near neighbor interactions and hopping parameters 
which break x-y symmetry. Although our results are 
for choices in which only hopping or interactions are 
anisotropic, we have also studied a few cases where 
x/y symmetry is broken in both terms^. These more 
complex situations also exhibited phase separation as 
opposed to supersolidity. Of course, we cannot exclude 
the possibility that a supersolid is stable in some other 
choice of parameter values. 

The motivation for studying such situations is to gain 
a more complete understanding of the interplay between 
real space particle arrangements and supcrflow in 
quantum systems. For strongly correlated fcrmions, this 
is a central puzzle of unconventional superconductivity: 
Do antiferromagnetic (AF) and charge order (stripes) 
compete or coexist with pairing? These issues arise not 
only in cuprate materials, which have AF ordering at 
q = (7r, 7r) but also in the more recently discovered iron- 
pnictides whose magnetic ordering is at q = (0, n): AF 
in one direction and ferromagnetic in the other. 

Since bosonic systems are more simple to study 
computationally than fermionic ones, it would be 
interesting to extend the present work to 'spin-1/2' 
bosons^i where one might have the analogs of both 
spin and charge order along with the possibility of 
superfluidity. 
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